# tested with python 3.2, numpy 1.6.1, scipy 0.10.1

import math
import csv;
from numpy import linspace
from scipy import exp, array
from scipy.signal import firwin
from scipy.fftpack import fft, ifft
import genblepcommon

class BiquadFilter(object):
    
    def __init__(self, b0, b1, b2, a1, a2):
        
        self.b0 = b0
        self.b1 = b1
        self.b2 = b2
        self.a1 = a1
        self.a2 = a2
        self.reset()

    def reset(self):

        self.x1 = 0.0
        self.x2 = 0.0
        self.y1 = 0.0
        self.y2 = 0.0

    def filter(self, x0):

        y0 = self.b0*x0 + self.b1*self.x1 + self.b2*self.x2 - self.a1*self.y1 - self.a2*self.y2
        self.x2 = self.x1
        self.x1 = x0

        self.y2 = self.y1
        self.y1 = y0

        return y0

def make_rc_lopass(sample_rate, freq):

    omega = 2 * math.pi * freq / sample_rate;
    term = 1 + 1/omega;
    return BiquadFilter(1/term, 0.0, 0.0, -1.0 + 1/term, 0.0);

def make_butterworth(fs, fc, res_db):

    # 2nd-order Butterworth s-domain coefficients are:
    #
    # b0 = 1.0  b1 = 0        b2 = 0
    # a0 = 1    a1 = sqrt(2)  a2 = 1
    #
    # by tweaking the a1 parameter, some resonance can be produced.

    res = 10.0 ** (-res_db / 10.0 / 2)
    ar = z_transform(1, 0, 0, 1, math.sqrt(2) * res, 1, fc, fs)
    return BiquadFilter(*ar)

# observe: a and b are reversed here. To be absolutely clear:
# a is the nominator and b is the denominator. :-/
def z_transform(a0, a1, a2, b0, b1, b2, fc, fs):
    
    # prewarp s-domain coefficients
    wp = 2.0 * fs * math.tan(math.pi * fc / fs)
    a2 /= wp * wp
    a1 /= wp
    b2 /= wp * wp
    b1 /= wp
    
    # compute bilinear transform and return it
    bd = 4 * b2 * fs * fs + 2 * b1 * fs + b0
    return [
        (4 * a2 * fs ** 2 + 2 * a1 * fs + a0) / bd,
        (2 * a0 - 8 * a2 * fs ** 2)           / bd,
        (4 * a2 * fs ** 2 - 2 * a1 * fs + a0) / bd,
        (2 * b0 - 8 * b2 * fs ** 2)           / bd,
        (4 * b2 * fs ** 2 - 2 * b1 * fs + b0) / bd,
    ]
           
def run_filter(flt, table):

    flt.reset()
    
    # initialize filter to stable state
    for i in range(10000):
        flt.filter(table[0])
    
    # now run the filter
    newtable = []
    for i in range(len(table)):
        newtable.append(flt.filter(table[i]))
    
    return newtable;
            
def main():
    
    a500_unfiltered, a1200_unfiltered = genblepcommon.prepare_unfiltered()

    # make digital models for the filters on Amiga 500 and 1200.
    filter_fixed5khz = make_rc_lopass(genblepcommon.AMIGA_CLOCK, 4900.0)

    # the leakage filter seems to reduce treble in both models a bit
    # the A500 filter seems to be well modelled only with a 4.9 kHz
    # filter although the component values would suggest 5 kHz filter.
    filter_leakage = make_rc_lopass(genblepcommon.AMIGA_CLOCK, 32000.0)

    filter_led = make_butterworth(genblepcommon.AMIGA_CLOCK, 3275.0, -0.70)
         
    amiga500_off = run_filter(filter_fixed5khz, a500_unfiltered)
    amiga1200_off = run_filter(filter_leakage, a1200_unfiltered)
    amiga500_on  = run_filter(filter_led, amiga500_off)
    amiga1200_on = run_filter(filter_led, amiga1200_off)

    genblepcommon.blep_calc('org-A500', amiga500_off, amiga500_on, a500_unfiltered);
    genblepcommon.blep_calc('org-A1200', amiga1200_off, amiga1200_on, a1200_unfiltered);

if __name__ == '__main__':

    genblepcommon.AMIGA_CLOCK = 3546895
    genblepcommon.VIDEO_MODE = 'PAL'

    main()

    genblepcommon.AMIGA_CLOCK = 3579545
    genblepcommon.VIDEO_MODE = 'NTSC'
    main()

